######测试贝叶斯代码
varis1<- c('MECPP_quantile','MnBP_quantile','MEHHP_quantile',
           'MEOHP_quantile','MiBP_quantile','cxMiNP_quantile','MCOP_quantile',
           'MCPP_quantile','MEP_quantile','MBzP_quantile')
demoB <-demo[,c('Age','Sex','race','edu','FPL','activity','smoke',
                'drink','Hypertension','CKD','urinecreatinine','BMI','URXECP', 'URXMBP', 
                'URXMHH', 'URXMOH','URXMIB','URXCNP','URXCOP', 'URXMC1','URXMEP','URXMZP',
                'Hyperuricemiaint','LBDSUASI')]
colnames(demoB) <- c('Age','Sex','race','edu','FPL','activity','smoke',
                     'drink','Hypertension','CKD','urinecreatinine','BMI','MECPP','MnBP','MEHHP',
                     'MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP','Hyperuricemiaint','LBDSUASI')

demoB$MECPP
demoTest1 <- demo
median(demoB$MECPP)
demoTest1$MnBP <- median(demoTest1$MnBP)
demoTest1$MECPP <- median(demoTest1$MECPP)
demoTest1$MEOHP <- median(demoTest1$MEOHP)
demoTest1$MiBP <- median(demoTest1$MiBP)
demoTest1$MiNP <- median(demoTest1$MiNP)
demoTest1$MCOP <- median(demoTest1$MCOP)
demoTest1$MCPP <- median(demoTest1$MCPP)
demoTest1$MEP <- median(demoTest1$MEP)
demoTest1$MBzP <- median(demoTest1$MBzP)


demoTest1 <- demoTest1[,c('Age','Sex','race','edu','FPL','activity','smoke',
                     'drink','Hypertension','CKD','urinecreatinine','BMI','URXECP', 'URXMBP', 
                     'URXMHH', 'URXMOH','URXMIB','URXCNP','URXCOP', 'URXMC1','URXMEP','URXMZP',
                     'Hyperuricemiaint','LBDSUASI')]
colnames(demoTest1) <- c('Age','Sex','race','edu','FPL','activity','smoke',
                     'drink','Hypertension','CKD','urinecreatinine','BMI','MECPP','MnBP','MEHHP',
                     'MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP','Hyperuricemiaint','LBDSUASI')
y <- demoTest1$Hyperuricemiaint
Z <- demoTest1[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP',
              'MCOP','MCPP','MEP','MBzP')]
X <- demoTest1[,c('Age','urinecreatinine')]
#X <- demoB[,c('Age','Sex','race','edu','FPL','activity','smoke',
#             'drink','Hypertension','CKD','urinecreatinine','BMI')]
fitkm <- kmbayes(y =y, Z = Z, X = X, iter = 2,verbose = FALSE, varsel = TRUE,family = 'binomial')
fitkm <- kmbayes(y, Z = Z, X = NULL, iter = 2, verbose = FALSE, varsel = TRUE, family="binomial",est.h=TRUE)
fitkm <- kmbayes(y,Z=Z,X=NULL,iter=2,verbose = FALSE,varsel = TRUE,family='binomial',est.h = TRUE)
pred.resp.univar <- PredictorResponseUnivar(fit = fitkm)
ggplot(pred.resp.univar, aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) + 
  geom_smooth(stat = "identity") + 
  facet_wrap(~ variable) +
  ylab("h(expos)") +
  xlab("Phthalate metabolites")




risks.overall = OverallRiskSummaries(fit=fitkm,qs=seq(0.25,0.75,by=0.05),q.fixed = 0.5)
risks.overall
ggplot(risks.overall,aes(quantile,est,ymin=est-1.96*sd,ymax=est+1.96*sd))+
  geom_hline(yintercept = 0,lty=2,col='brown')+
  geom_pointrange()









demo$MECPP_quantile <- cut(demo$URXECP,breaks = quantile(demo$URXECP),labels = c(1, 2, 3, 4))
demo$MnBP_quantile <- cut(demo$URXMBP,breaks = quantile(demo$URXMBP),labels = c(1, 2, 3, 4))  
demo$MEHHP_quantile <- cut(demo$URXMHH,breaks = quantile(demo$URXMHH),labels = c(1, 2, 3, 4))  
demo$MEOHP_quantile <- cut(demo$URXMOH,breaks = quantile(demo$URXMOH),labels = c(1, 2, 3, 4))  
demo$MiBP_quantile <- cut(demo$URXMIB,breaks = quantile(demo$URXMIB),labels = c(1, 2, 3, 4))  
demo$cxMiNP_quantile <- cut(demo$URXCNP,breaks = quantile(demo$URXCNP),labels = c(1, 2, 3, 4))  
demo$MCOP_quantile <- cut(demo$URXCOP,breaks = quantile(demo$URXCOP),labels = c(1, 2, 3, 4))  
demo$MCPP_quantile <- cut(demo$URXMC1,breaks = quantile(demo$URXMC1),labels = c(1, 2, 3, 4))  
demo$MEP_quantile <- cut(demo$URXMEP,breaks = quantile(demo$URXMEP),labels = c(1, 2, 3, 4))  
demo$MBzP_quantile <- cut(demo$URXMZP,breaks = quantile(demo$URXMZP),labels = c(1, 2, 3, 4)) 
demoTest1 <- demo
demoTest1 <- demoTest1[c(1:200),]
demoTest1 <- demoTest1[,c('Age','Sex','race','edu','FPL','activity','smoke',
                          'drink','Hypertension','CKD','urinecreatinine','BMI','URXECP', 'URXMBP', 
                          'URXMHH', 'URXMOH','URXMIB','URXCNP','URXCOP', 'URXMC1','URXMEP','URXMZP',
                          'Hyperuricemiaint','LBDSUASI','MECPP_quantile','MnBP_quantile','MEHHP_quantile',
                          'MEOHP_quantile','MiBP_quantile','cxMiNP_quantile','MCOP_quantile',
                          'MCPP_quantile','MEP_quantile','MBzP_quantile')]
colnames(demoTest1) <- c('Age','Sex','race','edu','FPL','activity','smoke',
                         'drink','Hypertension','CKD','urinecreatinine','BMI','MECPP','MnBP','MEHHP',
                         'MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP','Hyperuricemiaint','LBDSUASI',
                         'MECPP_quantile','MnBP_quantile','MEHHP_quantile',
                         'MEOHP_quantile','MiBP_quantile','cxMiNP_quantile','MCOP_quantile',
                         'MCPP_quantile','MEP_quantile','MBzP_quantile')
demoTest1 <- na.omit(demoTest1)
getwd()
library(bkmr)
library(ggplot2)
write.csv(demoTest1,'bkmrdemo.csv')
demoTest1 <- bkmrdemo
y <- demoTest1$Hyperuricemiaint
Z <- demoTest1[,c('MECPP_quantile','MnBP_quantile','MEHHP_quantile',
                  'MEOHP_quantile','MiBP_quantile','cxMiNP_quantile','MCOP_quantile',
                  'MCPP_quantile','MEP_quantile','MBzP_quantile')]
Z <- demoTest1[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP',
                  'MCOP','MCPP','MEP','MBzP')]
X <- demoTest1[,c('Age','urinecreatinine')]
for(i in 1:10){
  Z[,i] <- log(Z[,i])
}
Z$MECPP_quantile <- as.numeric(Z$MECPP_quantile)
Z$MnBP_quantile <- as.numeric(Z$MnBP_quantile)
Z$MEHHP_quantile <- as.numeric(Z$MEHHP_quantile)
Z$MEOHP_quantile <- as.numeric(Z$MEOHP_quantile)
Z$MiBP_quantile <- as.numeric(Z$MiBP_quantile)
Z$cxMiNP_quantile <- as.numeric(Z$cxMiNP_quantile)
Z$MCOP_quantile <- as.numeric(Z$MCOP_quantile)
Z$MCPP_quantile <- as.numeric(Z$MCPP_quantile)
Z$MEP_quantile <- as.numeric(Z$MEP_quantile)
Z$MBzP_quantile <- as.numeric(Z$MBzP_quantile)
set.seed(1000)
fitkm <- kmbayes(y, Z = Z, X = NULL, iter = 1000, verbose = FALSE, varsel = TRUE,family = 'binomial',est.h = TRUE)
# str(Z)
risks.overall = OverallRiskSummaries(fit=fitkm,qs=seq(0.25,0.75,by=0.05),q.fixed = 0.5)
risks.overall
ggplot(risks.overall,aes(quantile,est,ymin=est-1.96*sd,ymax=est+1.96*sd))+
  geom_hline(yintercept = 0,lty=2,col='brown')+
  geom_pointrange()


pred.resp.univar <- PredictorResponseUnivar(fit = fitkm)
ggplot(pred.resp.univar, aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) + 
  geom_smooth(stat = "identity") + 
  facet_wrap(~ variable) +
  ylab("h(expos)") +
  xlab("Phthalate metabolites") +
  xlim(-2,6)












